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Abstract 

Precise characterization of the mutual information of MIMO systems is required to assess the 

throughput of wireless communication channels in the presence of Rician fading and spatial correla- 
te 

tion. Here, we present an asymptotic approach allowing to approximate the distribution of the mutual 
information as a Gaussian distribution in order to provide both the average achievable rate and the 
outage probability. More precisely, the mean and variance of the mutual information of the separately- 
correlated Rician fading MIMO channel are derived when the number of transmit and receive antennas 
grows asymptotically large and their ratio approaches a finite constant. The derivation is based on the 
replica method, an asymptotic technique widely used in theoretical physics and, more recently, in the 
performance analysis of communication (CDMA and MIMO) systems. The replica method allows to 
analyze very difficult system cases in a comparatively simple way though some authors pointed out that 



X 



its assumptions are not always rigorous. Being aware of this, we underline the key assumptions made in 
this setting, quite similar to the assumptions made in the technical literature using the replica method in 
their asymptotic analyses. As far as concerns the convergence of the mutual information to the Gaussian 
distribution, it is shown that it holds under some mild technical conditions, which are tantamount to 
assuming that the spatial correlation structure has no asymptotically dominant eigenmodes. The accuracy 
of the asymptotic approach is assessed by providing a sizeable number of numerical results. It is shown 
that the approximation is very accurate in a wide variety of system settings even when the number of 
transmit and receive antennas is as small as a few units. 



* This work is supported by the STREP project No. IST-026905 (MASCOT) within the sixth framework programme of 
the European Commission. These results have been presented in part at the IEEE GLOBECOM 2006, Communication Theory 
Symposium. 
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I. Introduction 

During the last decade, multiple-input multiple-output (MIMO) wireless systems attracted 
much interest because of the prediction of outstanding capacity gains with respect to corre- 
sponding single-input single-output (SISO) systems [11], [46], [50]. 

Early work in this area was mostly based on the assumption of independent and identical 
distributed (iid) Rayleigh fading paths justified by the rich scattering assumption [11]. 

However, experimental and theoretical results showed that more realistic MIMO channel 
models are required to obtain more accurate results, which can dramatically impact the potential 
capacity gains [12], [13], [42]. Thus, more realistic channel models have been proposed in the 
literature accounting for both multipath correlation and the presence of a line-of-sight (LOS) 
component (Rician fading). 

In this work we focus on the separately-correlated MIMO Rician fading channel whose channel 
matrix can be written as 

H = H + R 1/2 H W T 1/2 , 

where H is a constant matrix accounting for the LOS component, R, T are the nonnega- 
tive definite receive and transmit correlation matrices, respectively, and H w is a matrix of iid 
circularly-symmetric zero-mean complex Gaussian random variables. We assume that the channel 
parameters (H, R, T) are unknown at the transmitter whereas the receiver knows the channel 
matrix realization exactly (i.e., following the terminology introduced by Goldsmith et al. [13], 
there is perfect channel-state information at the receiver (CSIR) and no channel distribution 
information at the transmitter (CDIT)). 

In the MIMO Rician fading channel literature, many works focus on the spatially uncorrelated 
channel. In particular, Hosli et al. [19] show that the mutual information is monotonically non- 
decreasing in the singular values of H. Jayaweera and Poor [22] derive the exact ergodic capacity 
when only the Rice factor K known at the transmitter by using the joint eigenvalue distribution 
of the noncentral Wishart matrix, and lower and upper bounds to the ergodic capacity when full 
CDIT is available. Hansen and Bolcskei [18] investigate the case of high SNR with unit-rank 
average channel matrix. Venkatesan et al. [47] derive the capacity-achieving covariance matrix 
of the channel, generalizing previous results relevant to the case of unit rank H [21]. Kang and 
Alouini [25] obtain the exact ergodic capacity by calculating the determinant of a matrix whose 
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entries contain confluent hypergeometric functions of the second kind and the exact capacity 
variance in terms of Meijers G-functions. 

Many recent works propose upper and lower bounds to the capacity of the separately cor- 
related Rician fading MIMO channel. McKay and Collings [31], [32] derive upper and lower 
bounds to the ergodic mutual information, based on the method of zonal polynomials, which are 
asymptotically tight as the SNR grows to infinity. They also provide an asymptotic (in the SNR) 
expression of the mutual information variance. The bounds have been further refined in [33]. In 
the same context, Cui et al. [3] obtain similar upper and lower bounds to the ergodic capacity 
by using a determinant expansion in terms of minors whose moments are calculated by use of 
hypergeometric matrix functions. In a recent paper, Lebrun et al. [28] use asymptotic bounds 
to confirm that the mean and variance of the MIMO channel capacity approach the values of 
the corresponding underlying scattering channel when the number of antennas is large and the 
average channel matrix has unit rank. 

Another way to circumvent the difficulties related to the derivation of an exact expression 
of the capacity consists of assuming an asymptotically large number of transmit and receive 
antennas. 

In this contest, Moustakas et al. [34] show that the capacity of the correlated Rayleigh 
fading MIMO channel converges to a Gaussian random variables whose mean and variance are 
calculated in closed form. Their derivation is based on approximating the moment generating 
function (MGF) of the capacity by using the replica method, a technique originally introduced 
in the context of statistical physics and successfully applied to several communications problems 
(see, e.g., [14], [37], [43], [49]). 

A considerable number of research works in this area are based on the use of Stieltjes 
Transform and lead to results compatible to those presented here. Among them, Dumont et 
al. [5] derived the asymptotic mutual information mean in the case of single-sided (receiver 
side) separately-correlated Rician fading. Hachem et al. [15], [16] derive the asymptotic mutual 
information mean in the case of uncorrelated (but not iid) Rician fading. Hachem et al. [17] 
derive the asymptotic mutual information mean and variance in the case of separately-correlated 
Rayleigh fading and consider the convergence of the random mutual information to the Gaussian 
distribution. 

Another approach to show that the asymptotic capacity distribution converges to the Gaussian 
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distribution is based on the study of the channel matrix singular values. Focusing on the corre- 
lated Rayleigh fading MIMO channel, Martin and Ottersten [30] show that, subject to certain 
conditions on the correlation between the channel elements, the limiting capacity distribution is 
Gaussian. 

Recent results in this area are aimed at extending the scenario to the MIMO communication 
channel with interference. This direction was taken originally by Moustakas et al. in [34] for 
the Rayleigh fading channel. More recently, Riegler and Taricco [40] derived the second-order 
statistics of the separately-correlated Rician fading MIMO channel mutual information in the 
presence of interference. This work was based on the use of Grassman variables and supermatrix 
theory (see references therein). 

A further extension of the results presented by this paper is the derivation of the ergodic 
capacity achieving covariance matrix. In this contest, the following works are worth mentioning. 

• Dumont et al. provide in [6] an asymptotic method (based on the Stieltjes transform) to 
derive the ergodic capacity achieving covariance matrix for the separately correlated Rician 
fading MIMO channel. The key ingredient is the asymptotic ergodic mutual information 
formula, which the authors refer to a preliminary version of [16], and is equivalent to 
the expression derived in this paper (by using a the replica method instead of Stieltjes 
transforms). The authors compare their results against those obtained by Vu and Paulraj, 
who used an interior point with barrier optimization method, and showed that the asymptotic 
method has a considerable advantage in terms of time efficiency. 

• Taricco and Riegler provide in [45] an algorithm for the derivation of the ergodic capacity 
achieving covariance matrix for the separately correlated Rician fading MIMO channel with 
interference. This work uses an asymptotic approximation of the ergodic mutual information 
for a separately-correlated Rician fading MIMO channel in the presence of interference. 

Finally, it is worth mentioning the recent result obtained by Riegler and Taricco [41] in the 
area of multiuser MIMO communications. This result consists in the derivation of the ergodic 
sum-rate achieving covariance matrices of a separately correlated Rician fading multiple access 
MIMO channel and of the resulting ergodic capacity region. The paper uses the replica method 
to derive the ergodic mutual information expressions required to upper bound the partial sum 
rates of the multiuser system. 

Summarizing, the main contributions of this paper are listed as follows. 
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• We consider an asymptotic setting where the number of transmit and receive antennas 
grows to infinity, while their ratio approaches a finite value. A closed-form expression of 
the asymptotic mean and variance of the mutual information of the separately-correlated 
MIMO Rician fading channel for arbitrary input signal covariance matrix is derived. These 
results depend on a pair of scalar parameters obtained by solving two polynomial equations. 
One of the main technical difficulties encountered, in addition to those addressed in [34], 
is in the need to deal with 2x2 block matrices in the saddlepoint approximation which do 
not split nicely as in the Rayleigh fading case [34]. 

• A detailed study of the cumulant generating function of the mutual information in the 
asymptotic regime leads to prove that it converges in distribution to a Gaussian random 
variable, provided that the correlation matrices T and R satisfy a mild technical condition. 
This condition is tantamount to requiring that the channel has no asymptotically dominant 
eigenmode. 

• From a theoretical point of view, the key issues of uniqueness of the fixed-point equation 
solution and of reality of the variance expression are addressed in two appendices. The latter 
result is important also because it allows to say that the replica-symmetric saddlepoint is 
stable. 

• To illustrate the analytic applications of the results obtained, the spatially uncorrelated Rician 
fading MIMO channel is considered and its ergodic capacity achieving covariance matrix 
is derived by using the asymptotic approximation. The results obtained are compared to the 
existing nonasymptotic literature. It can be noticed that the eigenvectors of the optimum 
covariance matrix can be obtained very easily from the average channel matrix (key result 
from [21]). Moreover, also the eigenvalues are derived — in the asymptotic setup — by 
using a fixed-point water-filling algorithm that is much simpler than the numerical techniques 
proposed in the literature to solve similar problems [24] (in the nonasymptotic case). 

• Finally, a set of numerical results is presented to validate the asymptotic analytic method. 
These results depart from a basic system scenario and change one of the following param- 
eters at a time: number of antennas, antenna ratio, spatial correlation (summarized in the 
base of the exponential correlation considered, common for the transmit and receive sides), 
Rice factor, and SNR. As far as concerns the ergodic capacity, it is shown that the relative 
error between the asymptotic and exact result (the latter derived by extensive Monte-Carlo 
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simulation) is always lower than 1%, even when the number of antennas is as low as one or 
two. The mutual information standard deviation asymptotic results are less accurate, with 
a relative error rising to a few percent units. However, using the asymptotic Gaussianity to 
derive an approximate outage mutual information, the standard deviation error is shown to 
have a negligible impact on this channel metric, yielding a relative error lower than 1% in 
the cases considered. 

A. Notation and basic results 

We denote (column-) vectors and matrices by lowercase and uppercase boldface characters, 
respectively. The imaginary unit is j = \f^\. The ath element of a vector x is (x) . The (a, fe)th 
element of a matrix A is (A) a6 . 5 ab = 1 if a = b and if a ^ b. The transpose of a matrix A 
is A T . The Hermitian transpose of a matrix A is A H . The Hermitian part of a square matrix A 
is defined as H(A) = i(A + A H ). The trace of a matrix is Tr(A) = ^2 a (A.)aa. The exponential 
of the trace of a matrix is etr(A) = exp(Tr(A)). The minimum and maximum eigenvalues of a 
matrix A are A min (A) and A max (A), respectively. The minimum and maximum singular values 
of a matrix A are <7 min (A) and <r max (A), respectively. The Frobenius norm of a matrix A is 
|| A ||; its square can be written as ||A|| 2 = Tr(AA H ). The Kronecker product of two matrices A 
and B is A <g> B; If A G C miXni and B G C m2X ™ 2 , then A <g> B G C mirn2Xri1 " 2 and we can write 
it in block matrix form as [(A) a6 B] for a = 1, . . . , mi and b = 1, . . . , n\. Among the properties 
of the Kronecker product we recall the following [20, p. 475]: i) A <g> (B <g> C) = (A <g> B) <g> C; 
H) lL(A fc ®B fe ) = (U k A fc )® (U k B fc ); Hi) (A®B) H = A H ®B H ; iv) (A^B)" 1 = A^B" 1 ; 
v) Tr(A®B) = Tr(A) Tr(B); vi) det(A®B) = det(A) n det(B) m if A G C mxm andB G C nxn . 
vec(A) is the column vector obtained by stacking the columns of A on top of each other from 
left to right. A 1 / 2 is the matrix square-root of the Hermitian nonnegative definite matrix A 
and is defined as UA 1 ^ 2 U H where A = UAU H is the unitary factorization of A [20, p. 414]. 
The notation x ~ A/" c (/x, R) defines a vector of complex jointly circularly-symmetric Gaussian 
random variables with mean value /x = E[x] and covariance matrix R = E[xx H ] — fifi H and its 
joint probability density function (pdf) is given by 

/(x) = det(vrR)- 1 exp[-(x - /x) H R- 1 (x - /u)] . 
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E[X] and Var[X] are the mean and the variance of the random variable X, respectively, (x) 



max{0,x}. Given the functions f(x) and g(x), we have, for x 
limsup^^ \f(x)/g(x)\ < oo and f(x) = Q(g(x)) if 



oo, f(x) = 0{g{x)) if 



< liminf 

x^oo 





< limsup 


/(*) 




x— >oo 





< oo. 



II. MIMO CHANNEL MODEL 

We consider a narrowband block fading channel with nj transmit and receive antennas 
characterized by the following equation: 

y = Hx + z. (l) 

Here, x G C™ TXl is the transmitted signal vector, H G C" rX " t is the channel matrix, z G C™ RXl 
is the additive noise vector, and y G C" rX1 is the received signal vector. 

We assume that the additive noise vector contains iid entries (z) a ~ A/" c (0, 1) for a = 1, . . . , n R . 

The channel matrix models separately (or Kronecker) correlated Rician fading so that it can 
be written as 

H = H + R 1/2 H W T 1/2 

where H represents the mean value and is related to the presence of a line-of- sight signal 
component in the multipath fading channel, the Hermitian nonnegative definite matrices T, R 
are the transmit and receive correlation matrices, and (H w ) ab ~ A/" c (0, 1) for a = 1, . . . , n R and 
b = 1, . . . , nj. The covariance between different entries of H is 

cov((H)jj, (H)i>j>) = ^ Pi'"- 



E[(R 1 /2) jfc (H lu ) fe ,(T 1/2 )„(R 1/2 )^,(Hj^(T 1 /^ 

k,e,k',£' 

^(R 1//2 )ife(T 1/2 )£j(R 1/2 )*/ fe (T 1/2 )^, 



k,( 

(R),:,;/(T 



30' ' 



A. Normalizations 

Assume that the input signal vector x has zero mearjj] and covariance matrix Q = E[xx H ]. To 
simplify notation, we define 

H 4 HQ 1 / 2 and T 4 Q l ' 2 TQ l '\ 

1 There is no point in having a nonzero mean input since the signal power would be greater but the mutual information would 
remain the same by the relationship J(x;y) = ft(x) — /i(x|y) = /i(x — fi^) — h(x — H x \y) where fi x — E[x] [2, Th.9.6.3]. 
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so that the input signal covariance matrix is implicitly accounted for into H and T. 
Then, the total received power is given byi 



E[||yf] = EfUHxfl + EdlR^H^T^xfl + Et 



z 



1 21 



= Tr[HQH H ] + Tr{E[H"RH M ,]T 1 ^ 2 QT 1 ^ 2 } + Tr[I„ R ] 

= ||H|| 2 + Tr(R)Tr(T) + n R . (2) 

According to eq. ©, the channel Rice factor (defined as the ratio of the received direct to 
diffuse power [9]) and the SNR are given by: 

IIHII 2 (K + 1) Tr(T) Tr(R) 

K = " " and p = - — — — — — . (3) 

Tr(R) Tr(T) n R 

In the special case of scalar input signal covariance matrix, Q = ql nj , we have: 

||H|| 2 n - _g n * t (a\ 

Tr(R) Tr(T) K + 1 Tr(R) Tr(T) W 



Remark II.l The two definitions © and © of the Rice factor are different when the input 
signal covariance matrix is not proportional to the identity matrix. Definition © complies with 
the common knowledge that the Rice factor represents the ratio between the line-of-sight and 
the scattered received signal power. It has the disadvantage of depending on the input signal 
covariance matrix, which may become an issue when the channel capacity is investigated. In 
that case, when a constraint on the total input signal power (Tr(Q)) is set, one has to resort to 
definition © even though it does not meet the common meaning of the Rice factor. 

B. Asymptotic setting 

Here we define the asymptotic setting assumed for the derivation of the mutual information 
mean and variance when nj, n R — ► oo and n T /n R — > k, (0 < k < oo). 

More precisely, we consider the sequence of parameters n R = 1,2,... and n T = \nn R ~\ and 
a corresponding sequence of deterministic matrices H, T, R, Q and random matrices H w of 
suitable dimensions. Additionally, 

• The Rice Factor K and the SNR p are constant as n R — > oo. 

2 We have: E[(H^RH w )y] = £ M E[(H W )^(R) W (H W )^] = Tr(R)5„ 
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• The correlation matrices are normalized by 

Tr(T) = n T and Tr(R) = n R . (5) 

• The matrices T and H are normalized by 

III. Channel mutual information and cumulant generating function 

It is well known [2] that the mutual information of a MIMO channel with channel matrix H 
and input signal covariance matrix Q is given by 

/(H) = lndet(I„ R + HQH H ) nats/complex dimension. 

The pdf of this random variable can be derived by its MGF, defined as 

G{y) = E[exp(-z//(H))] = E[det(I nR + HQH H ) " . 

Our goal is to derive an expression of the cumulant generating function (CGF) of /(H), defined 
as g(v) = mG(z/)E This allows us to derive the mean and the variance of /(H) as follows: 

E[/(x; y)] = -g'(0 + ) and Var[/(x; y)] = g"(0+). 

In the following, we resort to the replica method to obtain an asymptotic expansion of g(u). 

A. The replica method 

The replica method was originally used in the study of spin glasses [7]. Later, it found 
application in other research areas, such as neural networks, coding, image processing, and 
communications [39]. 

Many research works describe the physical meaning of the replica method (see, e.g., [14], 
[37], [39], [43], [49]). Here, we confine ourselves to a short review of its basic assumptions. 

According to [39], the replica method applies to a sequence of random variables Zn (typically 
representing the partition function of a physical system), converging in distribution to some Z^. 
The method is aimed at obtaining the so-called free energy of the system: 

T= lim - — E\\nZ N } = lim -—[ lim ^-BZKA \. (6) 

3 The cumulant generating function of a real Gaussian random variable with mean /i and variance a 2 is given by g(v) = 
— v[i + \v 2 a 2 . 
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The method is convenient when direct calculation is impossible (or very hard) whereas it is 
relatively easy to calculate the limit linijv^oo E[Z^] for positive integer v. In typical applications, 
the partition function can be expressed as a conditional average: Z N = E[Zfu; xjv)|xjv]. Here, u 
represents the system microstate and xjy is a set of N independent quenched^ random parameters. 
The method is based on the definition of Z a ^ N = E[Z(u a ; xjy)|xjv] as replicas of Z N , obtained 
by averaging with respect to the independent microstates u a for a = 1, . . . , v, conditionally on 
the quenched parameters. Hence, we can write: 



G(v, N) = E[Z^j 



a=l 



The validity of the replica method is subject to verification of the following assumptions: 

1) (Extension from positive integers). The limit G(v) = Hin^^oo G(v, N) is a smooth function 
of v in v — 0. This function is derived for positive integer v and extended to a right 
neighborhood of v — 0. 

2) {Interchange of limits). The limits in © can be exchanged, so that 

T = lim^-( lim -±-G(v,N)\. 

3) (Replica symmetry). The derivation of the limit 

Y\ Z(u a ; xat) 



lim E 



L a=l 



is based on the saddlepoint approximation (see Appendix |E] for a brief summary of this 
topic and [1], [23] for a deeper account) and on the symmetry of the stationary saddlepoint 
of ria=i l^( u ai x iv)l with respect to the replicated microstate arguments u a . 
In our application, the role of the partition function is played by 

Z nR = exp[-/(H)] = det(I nR + HQH")- 1 . 

Our goal is to determine the asymptotic series expansion in v of 

^)^lnE[Z n g=lnE[exp(-/(H))], 

4 A physical system is said to be in quenched disorder when some random parameters characterizing its behavior do not 
evolve in time and are then said quenched or frozen. Spin glasses are a typical example. Quenched disorder is in contrast to 
annealed disorder where all random parameters evolve in time [10], [39]. 
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and to show that only the first and second-degree coefficients survive in the limit. This is 
tantamount to saying that the mutual information converges in distribution to a Gaussian random 
variable. 

B. Derivation of the cumulant generating function 

Setting A = l v , C = I nR + HQH H , B = D = in eq. (O of Appendix E we can write: 

~g{y) 4 m{E[det(I„ R + HQHV]} 

= m( / E[etr{-7rU H (I nR + HQH H )U}]du). 

Next, by setting A = I„, C = I„ T , B = -D = Q 1/2 H H U in eq. ([19]) of Appendix E we obtain: 

etr(-7rU H HQH H U) = / etr[-7r(V H V + U H HQ 1/2 V - V H Q 1/2 H H U)]dV. 

JC"T X " 

Hence, we can rewrite the CGF as: 

~ g {y) = In ( f d\J [ E[etr{-7r(U H U + V H V + U H HQ 1/2 V - V H Q 1/2 H H U)}]rfvl. 
Now, we calculate the expectation by observing that 

Tr(U H HQ 1 / 2 V) = Tr(U H HV) + Tr[(T 1/2 VU H R 1/2 )H w ]. 

We obtain: 

g{v) = ln{ / dU I etr{-vr(U H U + V H V + U H HV - V H H H U) 

-7r 2 U H RUV H TV}rfv| 
after using the result (obtained by applying, e.g., eq. (fT9l of Appendix |B]): 

E[etr(AH u , - H^A H )] = exp(- 1| A|| 2 ). 
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Finally, by setting A = 7rU H RU and B = 7rV H TV in eq. ([23]) of Appendix E we obtain: 
g(y) = In ( / dU I dV lim / etr(WZ - aWW T + /?ZZ T ) 

etr[-vr(U H U + V H V + U H HV - V H H H U + WU H RU + ZV H TV)]^(W, Z)| 



In j f etr(WZ)d/i(W, Z) f du [ rfv 



cxp 



-vr(u H ,v^ 



lun R + W ® R 



In <^ / etr(WZ)def 



® H H I unj + Z ® T 
_! / L„ R + W ® R I„ ® H 




d/x(W, Z) 



(7) 



-I y ® H H I WT + Z ® T 

Here, the exchange of the limit for a, /3 j and the integral in U, V is allowed from the 
Dominated Convergence Theorem [27, Th. 5.30] since the absolute value of the first integrand, 
assuming W = W+j W and Z = Z +j Z for real W, W , Z, Z determined by the integration 
contour chosen, is dominated by 



/«exp 



vr(||U|| 2 + IIVI 



a 



W + 



7rU H RU - Zr 



2a 



7rV H TV + Wr 



2a 



which is plainly integrable over (U, V, W, Z) G C nRX " x C nTX " x M t/X,/ x R uxu . Similarly, the 
exchange of integration order is allowed by Fubini's Theorem [27, Th. 5.47] because the first 
integrand in C7]) is measurable as its absolute value is integrable [27, Prop. 5.53(c)]. Finally, 
after the change of integration order and the application of the limit, we used Property IA.1I and 
eq. (O of Appendix El 



C. Saddlepoint approximation 

In order to derive the asymptotic limit of © as n T , n R — ► oo, we use the method of saddlepoint 
approximation (see Appendix |E] for a summary of the main results and [1], [23] for a deeper 
account). 

To this purpose, we look for a stationary point of the integrand's logarithm in ©, namely, 



:W.Z) A Tri WZ) - In del 1 ' 



lvn„ + W (8) R 



-I„ ® H h 



l unT + Z ® T 
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To this end, we resort to the following expansion: 



00 ( — ]\ k - 1 
5[lndet(X)] = TrKX-^X)* 1 



k=l 



Then, assuming replica symmetry [34], we look for a stationary point of </>(W, Z) of the form 
W = wl u , Z = z\ v for some positive w, z. We have: 



(iuIi,, zl v ) — v < wz — In det 



I nR + wR 



H 



H I„ T + zT, 

The values of iu and z are obtained by setting the first-order variation of 0(W, Z) equal to zero. 
More generally, we can write the total variation at (wl u , zl v ) as: 

5(j)(wl v , zl u ) = Tr(w5Z + z5W + 5W5Z) 



00 / i\fc 



k=l 



k 



5W®(AiR) 5Z®(BiT 
5W®(dR) 5Z(8)(DiT 




(8) 



where [20, p. 18]: 



Ai 
Bi 



[I^+^R + H^ + ^T)- 1 ^]- 1 

(I„ T + zf )" 1 H H A 1 = -B" 

[I nT + + H H (I nR + wR)- 1 ^]- 1 



(9) 



Next, we focus on the second-order expansion of 5<p(wl L ,, z\ v ), since we show (Theorem lIV.ll) 
that subsequent terms vanish as np> — > 00. We have: 

^(tul,,, = [2 - Tr(AiR)] Tr(5W) + [w- Tr(DiT)] Tr(SZ) 

+i Tr[(AiR) 2 ] Tr(5W 2 ) + ~ Tr[(DiT) 2 ] Tr(5Z 2 ) 

+ [1 + Tr(BiTCiR)] Tr(5W<5Z). 

Thus, the stationary point is characterized by the values of w, z that are solutions of the 
following equations: 

w = Tr(DiT) =Tr|[^I nT + f- 1 + T- 1 H H (I nR +wR)- 1 H]- 1 
z = Tr(AiR) = Tr |[ioI nR + R 1 + R _1 H(I nT + 2;T) _1 H H ] _1 
The uniqueness of the solution of (fTOl) is proved in Appendix [Ql In the following, the solution 
of (flOl) will be denoted by (w,z), and the matrices defined in © after setting w = w,z = z 
will be denoted by Ai, Bi, Ci, Di, respectively. 



(10) 
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By the previous results, the saddle-point asymptotic approximation of g(y) (calculated along 
the directions of steepest descent from the stationary point) can be written as: 

I nR + wR H 



9\y) ~ v 



+ ln 




etr 



: |aX 2 -j3Y 2 + 2j 7 XY| dXdY^ 



where a = Tr[(AiR) 2 ], $ = Tr[(DiT) 2 ], and 7 = 1 + Tr(B 1 TC 1 R). Using the trace expansion 
Tr[d(X) 2 - /3(Y) 2 + 2 j 7 XY] = ]>>(X) 2 a - /3(Y) 2 a + 2 j 7 (X) aa (Y) aa ] 

a 

+2 ^[d(X) ab (X) ba - $(Y) ab (Y) ba + 2 j 7 (X) a6 (Y) 6o ], 

the CGF can be readily evaluated as 

I nR + w~R 



a<b 



g\y) ~ v 



wz 



— In det 



H H 



H 



It can be shown that < 7 2 — a/3 < 1 (see Appendix 10). As noticed in [36, Sec. IV-C] in 
a different setting, this inequality is required to guarantee the local stability of the replica- 
symmetric saddlepoint against variations around it. 

Thus, we have the following asymptotic expressions of the mean and variance of J(x;y): 



E[/(x;y)] ~ /i/ = lndet 



I nR + wK 
-H H 



H 



WZ 



(11) 



^ Var[/(x;y)] ~ aj 4 - ln( 7 2 - a$) 
expressed in (nat/complex dimension) and (nat/complex dimension) 2 , respectively. 

IV. Asymptotic Gaussianity 

The asymptotic Gaussianity of the mutual information requires that all terms of order k > 2 in 
the series expansion ® vanish as n R — > 00. This is shown under some mild technical conditions 
expressed in the statement of the following theorem. 



Theorem IV.l In the asymptotic setting specified in Section \II-B\ assuming further that the 
matrices R and T = T 1 / 2 QT 1 / 2 satisfy the conditions 



(TmaxlR) 



6(1) 



and 



ow(T) = eK- 1 ), 
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as — > oo, then the mutual information /(x; y) converges in distribution to a Gaussian random 
variable with mean and variance specified in (flTT) . 



Remark IV.l The conditions in Theorem IIV. 1 1 are tantamount to saying that the spatial cor- 
relation structure has no asymptotically dominant eigenmodes, i.e., both R and T do not 
have asymptotically dominant eigenvalues. Otherwise, if R and T had asymptotically dominant 
eigenvalues, the asymptotic nature of the MIMO channel would be drastically changed since only 
a small subset of antennas would play an asymptotically dominant role while the remaining ones 
would be asymptotically irrelevant. 

Considering T instead of T implies that also the structure of the input covariance matrix 
might affect the asymptotic evolution of the MIMO system. As an example, restricting power 
allocation to a finite subset of transmit antennas would be equivalent to nullify the corresponding 
columns of the channel matrix H and then would drastically change the asymptotic behavior of 
the MIMO channel. 

Remark IV.2 From the property of asymptotic Gaussianity, the outage probability corresponding 
to a given rate R can be approximated by 



with e denoting the outage probability. For example, I w % ~ [ii — I.2807 and I\% ~ /i/ — 2.33a/. 
As a result, it would be desirable to have a large mean fj,j and small standard deviation 07, the 
latter being a mere consequence of the fact that the outage mutual information improves as the 
cumulative mutual information distribution transition from to 1 becomes steeper. 

5 Here, Q(x) = P(7V(0, 1) > x) = -4^ exp(-u 2 /2)du. 
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Proof: See Appendix O 



P « P(M( m , aj)<R) = Qi Si 

as iV — ► oou Then, the outage mutual information is given asymptotically by 




(12) 
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V. Analytic example: Ergodic capacity of the uncorrelated channel 

Here we present an application of the asymptotic analytic results obtained to the calculation 
of the ergodic capacity of the spatially uncorrelated Rician fading MIMO channel. Assuming 
R = I nR and T = I„ T , the fixed-point eqs. (flOl) become: 

w = Tr|[zI„ T + Q- 1 + H H H/(l + ^)]- 1 } 
z = TrjKl + ^ + H^ + Q- 1 )-^"]- 1 } 
After some linear algebra, we can write the asymptotic ergodic mutual information as follows: 



fJ'i = n<R ln(l + w) + In det 



wz (14) 



where derives from H H H = U^A^U^ and Q = U^QU^. 

From Hadamard's inequality [20, Th. 7.8.1], we know that the maximum determinant of a 
positive definite diagonal matrix is upper bounded by the product of its diagonal elements and 
equality holds if and only if the matrix is diagonal. Thus, under the power constraint Tr(Q) = 
Tr(Q) < p, we can apply a standard water-filling argument [2] and obtain: 

V 1 +w + {A H )uJ + 
where (x) + = max{0, x} and £ can be obtained by solving the equation: 

l + w + (A H ) u J 



Notice that eqs. (1131) and (1161) have to be solved simultaneously since we have a mutual in- 
terdependence between Q and the pair (w,z). In other words, the solution can be derived by 
implementing a simple iterative fixed-point water-filling algorithm. 

Finally, we notice that the structure of the asymptotic ergodic capacity achieving covariance 
matrix is consistent with [21], which showed that the capacity achieving covariance matrix and 
the matrix H H H have the same eigenvectors, whereas the problem of calculating the eigenvalues 
requires numerical optimization techniques such as those used in [24]. It is then clear the 
advantage of the asymptotic approximation, which allows to derive analytic results very simply. 
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VI. Numerical examples 

In this section we consider a baseline system scenario where we assume that the average 
channel matrix is given by (H)y = 1 for all i, j, the spatial correlation matrices are of exponential 
type with common base a, namely, (T)^ = (R)ij = a'* - - 7 ', there is no CDIT, so that capacity 
is achieved by setting Q = ql t , and the specification is completed by the following set of 
parameters: 

' n R =4 
nj/n R = 1 

• (IV) 

K = 10 dB 

a =0 

Then, we consider the impact of changing each parameter in turn, with the aim of illustrating 
the accuracy of the asymptotic approximation proposed. For each parameter we plot the mean 
and standard deviation of the mutual information by using the asymptotic method (solid lines) 
and Monte-Carlo simulation (markers). 

A. Impact of the number of antennas 

Figures \T\ and |2] describe the impact of the number of antennas and of the SNR on the 
asymptotic approximation accuracy. It can be noticed that the approximation of the mean (er- 
godic mutual information) is always excellent (the maximum relative error being less than 2%, 
corresponding to the case of nj = np> = 2 and SNR = 30 dB). The accuracy of the standard 
deviation is very good when the SNR is or 10 dB (relative error always less than 1%) but only 
fairly accurate when the SNR is 20 or 30 dB (maximum relative error around 5%). However, it 
must be noted that the impact on the outage mutual information approximation (fl"2l) is attenuated 
by the fact that the standard deviation is considerably smaller than the mean when the SNR is 
sufficiently large. 

B. Impact of the antenna ratio 

Figures [3] and 0] describe the impact of the TX to RX antenna ratio and of the SNR on the 
asymptotic approximation accuracy. It is assumed that the number of receive antenna is fixed, 
n R = 4, and the number of transmit antennas ranges from nj — 1 to 10. Both the mean and 
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the standard deviation of the mutual information display very good accuracy with a maximum 
relative error smaller than 2%. 

C. Impact of the Rice factor 

Figures [5] and [6] describe the impact of the Rice factor and of the SNR on the asymptotic 
approximation accuracy. It can be noticed that the ergodic capacity accuracy is always very 
good in the cases observed (maximum relative error around 1%) while the standard deviation is 
slightly overestimated for low K (maximum relative error around 4%) and underestimated higher 
K (the threshold depending on the SNR). This effect can be explained by considering that the 
low K condition entails a larger amount of randomness in the channel matrix that increases the 
asymptotic variance expression. The deviation is larger for high SNR, as in the previous cases 
observed. Again, the error on the standard deviation has a modest relative impact on the outage 
mutual information which is the goal of the variance approximation. 

D. Impact of spatial correlation 

Figures [7] and [8] describe the impact of spatial correlation of the SNR on the asymptotic 
approximation accuracy. The relative error remains below 1% for the ergodic capacity and 3% 
for the standard deviation. The relative error increases as the SNR increases and correlation 
decreases. 

It is interesting to note that spatial correlation has a minor effect at moderate SNR, whereas 
it becomes important as the SNR increases. 

Summarizing the numerical results obtained so far we can say that randomness tends to 
increase the relative error between the analytic asymptotic approximation and the actual value 
(obtained via Monte-Carlo simulation). So, increasing the SNR and decreasing the correlation 
or the Rice factor produces an increase of randomness that entails a larger relative error. Besides 
this, it can be noticed that the relative error gets larger as the antenna ratio gets closer to 1. 
Finally, the relative error decreases as the number of antennas increases but this is a trivial 
consequence of the fact that the approximation is asymptotic in the number of antennas. 
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E. Outage mutual information 

In order to assess the accuracy of the proposed asymptotic approximation we consider here 
the outage mutual information, defined implicitly as follows: 



We compare the asymptotic analytic approximation (fl"2l) against the exact numerical result 
obtained by accurate (though lengthy) Monte-Carlo simulation. The results are illustrated in 
Figs. [9] and [TOl reporting I w % (i.e., the 10% outage mutual information) versus the number 
of antennas and the exponential correlation base a, respectively. The maximum relative error 
is around 1% in the cases considered, with the maximum corresponding to the case of 2 x 2 
MIMO with the largest SNR. This is comparable with the relative error already found for the 
ergodic capacity approximation. 



This paper presented an analytic approach, based on the replica method, allowing to approx- 
imate the statistics of the mutual information for a separately correlated Rician fading MIMO 
channel. We showed that the mutual information statistics approach the Gaussian distribution as 
the number of transmit and receive antennas grows large, with their ratio approaching a finite 
constant. 

More specifically, we saw that the mutual information mean (corresponding to the ergodic 
mutual information of the channel) yields a very accurate approximation of the real value (which 
was obtained by extensive Monte-Carlo simulation) with a relative error never larger than a few 
percent units. This remarkable accuracy was obtained not only with a large number of antennas 
but also in cases when the number of antennas is definitely small, even in the limiting SISO case! 
As a result, the asymptotic analysis becomes a valuable tool to assess the system performance 
whenever this has to be done analytically, such as in the case of covariance optimization with the 
goal of finding the channel capacity. This application is illustrated in the paper by considering 
the uncorrelated Rician fading channel. 

We also showed that the mutual information standard deviation approximation is very good. 
This is important in view of the application of these results to the derivation of the outage mutual 
information. We showed that the outage mutual information, obtained by using the asymptotic 




R>0 



(18) 



VII. Conclusions 
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mean and standard deviation, is very accurate in the cases considered, since the relative error 
is always lower than 1%. This result depends on the asymptotic Gaussianity of the mutual 
information, which is proven in the paper. 

Summarizing, the analytic approach presented allows to address concisely the performance of 
realistic models of MIMO channels (based on separately correlated Rician fading). The accuracy 
obtained is sufficient for most applications and eliminates the need to run computer intensive 
simulations or implement numerical methods in orderto optimize the system performance. Many 
applications of this asymptotic method are appearing in the conference literature. Among them, 
it is worth mentioning the optimization of the transmitted signal covariance, with and without 
interference, and the derivation of the capacity region of the multiuser multiple access MIMO 
channel. A simple example of the former (relevant to the case of uncorrelated Rician fading 
without interference) is given in the paper. 
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Appendix 

A. Property of the Kronecker product 

Property A.l For any matrices A, X, B,Y such that the product AX H BY exists and is a 
square matrix, we have: 

Tr(AX H BY) = vec(X) H (A T ® B)vec(Y). 
Proof: Assume that X = (xi, . . . , x m ) and Y = (yi, . . . , y n ). Then, 



Tr(AX H BY) = Tr 



t xf B yi ■ ■ • xf By„ \ 



A 



V 



x^B Yl ••■ x^By n y_ 
^(A) a , fc xfBy a = (A) a , 6 By a 

a,b b a 

vec(X) H (A T ®B)vec(Y). 



B. Integrals 

Integral B.l For every pair of Hermitian positive definite matrices A e c mxm ; q ^ £_nxn^ an( ^ 
for any matrices B, D G C mxn , we have: 

[ etr[-7r(AU H CU + B H U + U H D)]rfU = det(A T ® C)" 1 etr(7rA- 1 BC" 1 D H ) 

= det(A)~ n det(C)- m etr(7rA- 1 BC~ 1 D H ). 

(19) 



Proof: From Property IA.ll we have: 

Tr(AU H CU) = vec(U) H (A T ® C)vec(U). 
Then, we can write the integral on the lhs of (fT9l) as 

/ exp[-7r(u H Pu + b H u + u H d)]du 

j£mn 

where P = A T ® C, b = vec(B), and d = vec(D). Since A and C are Hermitian and positive 
definite by assumption, so are A T and P = A T ® C. 
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Applying the change of variables v = P 1 / 2 u, we obtain: 

(det P)- 1 / exp[-7r(v H v + b H P~ 1/2 v + v H P~ 1/2 d)]dv. 
This integral can also be written as 

ran „ 

(detP)" 1 ]! / exp[-n(\v a \ 2 + ~b* a v a + v*J a )}dv a 

a=l Jc 

where we let b = P x / 2 b and d = P x / 2 d. Setting v a = v ax + jv ay , we have 

\v a \ 2 + Kv a + v*d a = v 2 ax + v ax {b* a + d a ) + v 2 ay + jv ay (b* a - d a ). 

By using repeatedly the standard Gaussian integral 

/ exp[— 7r(x 2 + ax)]dx = exp(7ra 2 /4), a G C, 
Jr 



we obtain: 



/ exp[-vr(v H v + b H P~ 1/2 v + v H P~ 1/2 b)]dv 

J£mn 

mn 

J] { ex ^& + J «)V 4 ] exp[-vr(fe: - d a ) 2 /4] 



a=l 
mn 



= Yl exp(vr6* J a ) = exp(vrb H d) = exp(7rb H P _1 d). 

a=l 

Finally, since P 1 = (A -1 ) 7 ® C~\ applying Property lA.ll we get: 

vec(B) H ((A- 1 ) T <g> C- 1 )vec(D) = ^(A^BC^D"), 
which in turn leads to (fl9l). 



Remark .1 Ref. [34, eq.(124)] reports the following incorrect result: 

[ dr exp(— rt>) / exp(rt — ut) — - = exp(— uv) u,v E C. (20) 
Jr J] m 27rj 

The previous result is wrong because the function exp(— rv) is not a test function for the 6- 

sequence 8a{t — u) = f \ exp(j (r — u)t)^ (unless v is purely imaginary) since it is unbounded 

over the integration domain of r [26]. However, if we smoothen the test function by it by a 

suitable Gaussian factor, we obtain: 

[,[ f 2 n 2 1 dt 1 /— uv — au 2 + (3v 2 \ 
/ dr exp < —ar + tit + rt — ut — rv > = = exp . 

(21) 
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This is a Lebesgue integral provided that a, (3 > and it holds for any r ,t e M. Therefore, 
we can replace eq. (l20l) by: 

exp(— uv) = lim / dr [ exp { — ar 2 + (3t 2 + rt — ut — rvX . 

a,/3|Oj R+jro Jt +]R 1 J 27TJ 

The previous result allows us to prove the following proposition. 



(22) 



Integral B.2 (Robust Hubbard-Stratonovich transformation) For every pair of matrices A e 
C mxn ,B G C nxm , we have: 

etr(-AB) = lim / etr(-«WW T + (3ZZ 1 + WZ - AW - BZ)d//(W, Z) (23) 

where d//(W, Z) = Yla=i YlT=i ^jd(W) ab d( y Z) ba and integration is carried out along contours 
parallel to the real axis for each (W) ab and to the imaginary axis for each {Z) ba . 



Proof: The proof stems from repeated application of (1221) . 



C. Proof of Theorem \IV,1\ 

Proof: Using the complete expansion in ([8]), we can write the CGF as 



g(v) ~ v 




I nR + wR 



H 



+ ln 



-H H l nj + zT 



^«U 2 -/3V 2 + 2j'7UV 



dudv y 



U®(AiR) jV®(B x T) 
U®(C X R) jV®(DxT) 

Thus, in order to prove the convergence of g(u) to the Gaussian CGF, we have to show that, 



as n R — > oo, 



etr 



!<jdU 2 - / 3V 2 + 2j'7UV 



Tr 



U®(AiR) jV®(BiT) 
U®(CiR) jV<8>(DiT) 



dVdV-^0 (24) 



Here we used the series expansion (l + x)^ 1 exp(x-x 2 /2) = l + ^fcls c k xk - 1 - §a: 3 + ;j£ 4 - §x 5 + fx 6 - ff 2 7 + ■ 
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for A; > 3. Now, in order to balance the asymptotic order of the matrix multipliers of U and 
V, it is convenient to apply the change of variables U = 72r _1 / 2 TJ and V = n R l l 2 ~V. Then, the 
limit (|24)) becomes: 



etr 



|<jn R - 1 dU 2 -n R /3V 2 + 2j'7UV 



■Tr 



TJ ® (n R - 1 / 2 A 1 R) jV ® (n R 1 / 2 B 1 T) 
U ® (n R - 1 / 2 C 1 R) j V <8> (n R 1 / 2 D 1 T) 



c/Uc/V 



(25) 



and since 



a 



P 



Tr[(A x R) 2 ] < Tr(R 2 ) = 0(n R ), 
Tr[(DiT) 2 ] < Tr(T 2 ) = O^ 1 ) 



by the theorem's assumptions, it is plain to see that condition (1251) is equivalent to the limit 



Tr 



riR— >oo 







(26) 



U ® (n R - 1 / 2 A 1 R) jV ® (n R 1 / 2 B 1 T) 
TJ ® (riR-^CiR) jV ® (nR^DiT) 

To prove (1261) . we investigate the asymptotic order of the singular values of the matrices A X R, 
BiT, CiR, and DiT. We shall use the following linear algebra inequalities [20]: 

cr max (MN) < o"max(M)<j max (N), for any M, N with compatible sizes 
cr max (M) < cr max (N), for any square M, N with M < N 

Tr(M) < Tr(N), for any square M, N with M < N 

From the asymptotic setting definitions of Section Hl-BI we have the following results. 
1) Since 

AiR = [wl nR + R 1 + R^HQ^ + zT)- 1 !!"]- 1 , 
we have A X R < R, so that 



(Tmax(AiR) < a 

max 

(R) = e(i; 



and 



z = Tr(AiR) < Tr(R) = Q(n R ). 



2) Since 



DiT = [zl nj + T- 1 + T^H^ + wR)- 1 ^]- 1 , 
we have DiT < T, so that 

Cmax(DiT) < <x max (T ) = e^R- 1 ) and w = Tr(DiT) < Tr(T) = 6(1). 
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3) By using the previous inequalities on the maximum singular values we obtain: 

Tr[(AiR) fc ] = 0{n R ) and Tr[(D 1 T) fe ] = O(V^) 

which confirms, for k = 2, that a = 0(n R ) and (3 = 0(n R ~ r ). 

4) Setting H = (I„ R + u?R)" 1/2 H, by the definition of Bi in ©, we have: 

a max (B!T) = a max [(I nR + iZ;R)- 1/2 H(I nT + zT + H H H)- 1 T] 

< ^max[(I„ R + tfR)- 1/ V ln ax[H(I nT + + H H H) _1 ] <7 max (T) . 

Since (plainly) cr max [(I nR + u)R) _1 / 2 ] = 0(1) and 

aLjH(In T + ^T + H H H)- 1 l = A max [H H H(I nT + iT + H H H)- 2 l 



Hf±\-2l 



< A max [H H H(I nT + H H H 

Ai(H H H) ^ 1 
= max — — — — < - , 

* [1 + Aj(H H H)] 2 ~ 4 

we have 

C r max(BiT) = O^r- 1 ). 

5) Setting H = H(I„ T + zT)' 1 / 2 , by the definition of Ci in ©, we have: 
^(CiR) = a max [(I nT + ^T)- 1 / 2 H H (I nR + ri;R + HH H )- 1 R] 

< <w[(I nT + IT)- 1 / 2 ]^^"^ + wK + HH H )-Vmax(R). 
Since (plainly) a max [(I„ R + zT)- 1 / 2 ] = 0(1) and 

^[H^I^ + wR + HH")- 1 ] = A max [HH H (I nR + u»R + HH H )~ 2 ] 

< A max [HH H (I nR +HH H )- 2 ] 

A t (HH H ) ^ 1 

= max < - , 

* [1 + Aj(HH H )] 2 ~~ 4 

we have 

^m(CiR) = 0(1). 

In order to prove (|26l) . notice that the maximum singular values of the matrices on the rhs of 
the Kronecker products in 

*U® (fiR-^AiR) j V ® (nR^BiT) 
U ® (nR-V^CiR) jV ® (n R 1 / 2 D 1 T) 
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are all 0(n R l l 2 ) as n R — ► oo. Thus, if we define 

W(U,V)®M n , fc 0i 2 , fc (U,V)<8>M 12 , fc \ ^ /u®(n R -V2 Al R) jV ® (n R 1/2 BiT) 
W(U,V)®M 21 , fc 22 ,A ; (U,V)®M 22ife y yu<8> (n R ~ 1/2 CiR) jV ® (V^T) 

we have o" max (Mj Jjfc ) = 0(n R ~ fc / 2 ). Therefore, the trace 

Tr(0u lfc (U, V)) Tr(Mn, fc ) + Tr(0 22 , fc (U, V)) Tr(M 22jfc ) = O^r 1 ^ 2 ) 

approaches as n R — >■ oo for all k > 3, which proves (|24)) and the fact that g{y) converges to 
a Gaussian CGF. 



D. Uniqueness of the solution of eqs. (1101) 

The proof reported here was inspired by a similar proof included in [17] 
Proof: Eqs. (flOl) can be written equivalently as follows: 

Tr(DiT) 



1 



1 



w 

Tr(AiT) 



Tr {[wzl nj + wT- 1 + ibT^H"^ + wR)- 1 !!]- 1 
Tr \[wz! nR + zK 1 + zR- x H(I nT + ^T)" 1 ^]" 1 



(27) 



First, we calculate the partial derivatives of these functions with respect to w and z: 

< 



dipt 


Tr{DiTD![I nT + zT + H H (I nR 


+ «jR)- 2 H]} 


dw 


_ w 2 




dtp! 


Tr{(DiT) 2 } ^ 




< dz 


w 




dip 2 


Tr{(AiR) 2 } ^ 




dw 






dip 2 


Tr{AfRA x [I nR + wK + H(I nT 


+ zf)- 2 H H ]} 


^ dz 


z 2 





(28) 



< 

Then, we notice that ^(O, z) = oo, t/>x(oa,z) = 0, and dip\/dw < 0. Hence, xj>i(w,z) is a 
continuous monotonically decreasing function in w and the equation ipi(w,z) = 1 has a single 
solution w = g(z) E [0, oo), which is a continuous function of z by the Implicit Function 
Theorem [27, Th. 3.16]. Moreover, from the first of (1271) and from the first two inequalities of 
(1281) . we obtain: 



dipi I dz 
dtpi I dw 



< 0. 
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Then, the uniqueness of the solution of (flOl) requires that the equation h(z) = $2(9(2), z) — 1 
has a single solution. This result stems from the fact that ^(^(O), 0) = +00, $ 2 (g(oo), 00) = 0, 
and 



h'(z) 



■0 0) + 



dw 



dz 



< 0. 



w=g(z) 



The last inequality is equivalent to 

dtpi dip2 dipi dip2 
dw dz dz dw± w=g{z) 

which can be checked by direct substitution of the relevant expressions. 



> 



E. Saddlepoint approximation 

Here we report some basic facts about saddlepoint approximation. A full account on the subject 
can be found in [23]. In its simplest form, saddlepoint approximation deals with integrals of the 
type 

I— g(x) exp[—\h(x)]dx 
Jr 

in the limit for A — > 00. Assuming g(x) bounded and h(x) with a global minimum at x = xq, 
the integral can be approximated, for A — ► 00, as: 



g(xo) exp[— Xh(xo) — \h"(xo)(x — xq) 2 /2}dx 



2tt 



Xh"(x ) 



g{x )exp[-\h(x )}. 



A more refined derivation yields 



exp[-Xh(x )] I , ^ + g(5h"' - 3h" 2 h iv ) - 12g'h"ti" + 12g"h" 2 } 



X = Xq 



y/Xh"(x )/(2n) ' 24/i"3 

This result can be extended to a multidimensional contour integral 

g(x) exp[— A/i(x)]dx. 

Assuming \g(x)\ bounded and /i(x) real, smooth, and with a global minimum at x = x , we 
have, for A — » 00: 

5f(x )exp[-A/i(x )] 



/ « J 5f(x )exp[-A/i(x ) - A(x-x ) T H (x-x )/2]rfx 



Vdet[AH /(27r)] ' 

where H is the Hessian matrix of h(x) at x = x , which is assumed to be positive definite. 
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F. Proof of the inequalities < — af3 < 1 

Proof: First, we show that < 7 < 1, so that, since a/3 > 0, we have 7 2 — a/3 < 1. 

• The inequality 7 < 1 derives from 

7 = 1 + Tr^T CiR) = 1 - T^R^C^TCxR 1 / 2 ) = 1 - HT^CxR 1 / 2 )! 2 . 

• The inequality 7 > derives from 

7 = 1 + Tr(BiTCiR) 

= 1 - TrfAiH^ + zT)- l T(l nj + ST^i&MR] 

(i) ~ ~ ~ 

> l-^TrfAiH^ + iTj-WAiR] 

( => l-i- 1 Tr[(I nR -A 1 (I rtR +u ; R))A 1 R] 
Tr(A 2 R) + ?M 

z 

from (i) the matrix inequality (I„ T + zT) _1 T < z _1 I nT ; (ii) the definition of A^ and (Hi) 
the definition of a — Tr(AiR). Now, in order to handle the case of singular R, we define 
the positive definite matrix R e = eI„ R + R for e > 0. We have: 

Tr(A 2 R e ) — eTr(A 2 ) + wa 

7 > : 

z 

W Wg£j Tr(AiR e AiR e ) - 6 Tr( A 2 ) + 

~~ z 

(«) A max (R) _1 + t« v 

— > a 

z 

> 0, (29) 

from (i) the inequality Tr(UV) > A m i n (U) Tr(V) (holding for nonnegative definite matrices 

U, V as a direct consequence of [29, 9.H.l.h]); and (ii) taking the limit for e J. of the 

previous lower bound and the limit lim e j_ A min (R e ) = A max (R) _1 . 

The second part of the proof can be accomplished by combining the inequality 

- . A max (R)^ 1 + w v w v 
7 > a > —a, 

deriving from ([291) , and 

7 > Amax(T) " 1 + i /3 > 
which can be obtained in a similar way. Multiplying the two inequalities yields -y 2 — a/3 > 0. ■ 
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Fig. 1. Capacity mean versus number of antennas for different SNR's. 
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Fig. 2. Capacity standard deviation versus number of antennas for different SNR's. 
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Fig. 3. Capacity mean versus the TX to RX antenna ratio for different SNR's. 
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Fig. 4. Capacity standard deviation versus the TX to RX antenna ratio for different SNR 's. 
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Fig. 5. Capacity mean versus the Rice factor K (expressed in dB) for different SNR's. 
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Fig. 6. Capacity standard deviation versus the Rice factor K (expressed in dB) for different SNR's. 
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Fig. 7. Capacity mean versus exponential spatial correlation base a for different SNR's. 
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Fig. 8. Capacity standard deviation versus exponential spatial correlation base a for different SNR's. 
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Fig. 9. Outage capacity versus number of antennas for different SNR's. 
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Fig. 10. Outage capacity versus exponential spatial correlation base a for different SNR's. 
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